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Electron energy loss spectroscopy (EELS) has emerged as a powerful tool for the investigation of 
plasmonic nanoparticles, but the interpretation of EELS results in in terms of optical quantities, 
such as the photonic local density of states, remains challenging. Recent work has demonstrated 
that under restrictive assumptions, including the applicability of the quasistatic approximation and 
a plasmonic response governed by a single mode, one can rephrase EELS as a tomography scheme for 
the reconstruction of plasmonic eigenmodes. In this paper we lift these restrictions by formulating 
EELS as an inverse problem, and show that the complete dyadic Green tensor can be reconstructed 
for plasmonic particles of arbitrary shape. The key steps underlying our approach are a generic 
singular value decomposition of the dyadic Green tensor and a compressed sensing optimization for 
the determination of the expansion coefficients. We demonstrate the applicability of our scheme for 
prototypical nanorod, bowtie, and cube geometries. 
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Introduction 

Electron energy loss spectroscopy (EELS) is a powerful 
tool for the investigation of plasmonic nanoparticles.^i^ 
EELS is a technique based on electron microscopy and 
measures the probability of a swift electron to lose part of 
its kinetic energy through plasmon excitation as a func¬ 
tion of electron bea m po sition. Eollowing first proof of 
principle experiments,^^ in the last couple of years EELS 
has been exhaustively used for the investigation of plas¬ 
mon modes in single and coupled nanoparticles. 

Despite its success, the interpretation of EELS data 
in terms of optical quantities, such as the photo nic lo¬ 
cal density of state^(LDOS), remains challenging.!^ To 
overcome this problem, in Ref. |8] we formulated EELS 
as a tomography schem^ and showed that under cer¬ 
tain assumptions a collection of EELS maps can be 
used to reconstruct the three-dimensional mode profile 
of plasmonic nanoparticles. A similar approach was pre¬ 
sented independently by Nicoletti and coworkers,^ who 
demonstrated the applicability of the scheme for a sil¬ 
ver nanocube. Extracting three-dimensional information 
throu gh s ample tilting was also shown for a split-ring res¬ 
onator^ and a nanocrescent using cathodoluminescence 
imaging.l^ 

The problem with EELS tomography is that the mea¬ 
surement signal (the loss probability) is not simply the in¬ 
tegral of local losses along the electron trajectory, but in¬ 
volves a two-step process where the swift electron first ex¬ 
cites a particle plasmon and then performs work against 
the induced particle plasmon field. This leads to a non¬ 
local response function, which allows for a tomographic 
reconstruction only under restrictive assumptions, such 
as the applicability of the quasistatic approximation or 
a plasmonic response governed by a single mode. In this 
paper we use additional pre-knowledge, namely that the 
particle plasmon fields are solutions of Maxwell’s equa¬ 
tions and that the dyadic Green tensor^ can be decom¬ 


posed into modes, in order to rephrase EELS in terms 
of an inverse problem. We develop a rather generic 
model for the EELS probabilities, which depends on a 
few parameters, and determine the parameters such that 
the model data match as closely as possible the mea¬ 
sured data. Within this approach we are able to obtain 
most accurate reconstructions of the dyadic Green tensor, 
which, in turn, allows us to extract the three-dimensional 
photonic LDOS from a collection of tilted EELS maps. 
We demonstrate the applicability of our scheme for pro¬ 
totypical nanorod, bowtie, and cube geometries. 

Theory 

We start by analysing EELS within a semi-classical 
framework,!^] where a swift electron propagating with ve¬ 
locity V loses a tiny part of its kinetic energy by perform¬ 
ing work against the electric field E[re{t)] produced by 
itself. Eor sufficiently large velocities we can ignore ve¬ 
locity changes in the electron trajectory re(t) « Rq + vt, 
with Rq being the impact parameter. It is convenient to 
split E = Elbuik + ^^surf into a bulk contributioiP^ Elbuik, 
corresponding to the electric field within an unbounded 
homogeneous medium, and a surface contribution E^uri, 
corresponding to field modifications (including surface 
plasmons) from the interfaces between different materi¬ 
als. Bulk losses are due to Cherenkov radiation and elec¬ 
tronic excitations,® and the loss probability is obtained 
by simply multiplying the loss probability per unit length 
Tbuik(‘^)’ inside material j and for loss energy Hu, with 
the path length £j of the electron inside material j, 

rbuik(w) = ■ (1) 

3 

Bulk losses can be interpreted in terms of local scatter¬ 
ings where the electron emits a photon or excites elec¬ 
trons in the dielectric material, and loses part of its ki- 
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netic energies. To compute the surface loss probability, 
we integrate the work dW = eEguif ■ vdt performed by 
the electron over the entire trajectory, and decompose it 
into the different loss energies huj according to 

/ OO pOO 

V ■ Esurl[re{t)]dt = / hujTsnTl{i^)duj . (2) 
-OO J 0 

Thus, the energy loss probability becomeJ^ 


= :^ j R,e|e ■ EsnH[re{t),uj]^ dt, 

( 3 ) 

where we have explicitly indicated the dependence on the 
electron propagation direction and the impact parame¬ 
ter through Ri, = {v,Ro). To understand the physical 
process underlying Eq. it is convenient to introduce 
the current distribution J{r,t) = —evS{r — Vf-it)) of the 
swift electron and the dyadic Green tensoJ^ G{r,r',uj) 
that relates for a given frequency uj a current source at 
position r' to an electric field at position r via E{r,u;) = 
iijj^oG{r,r',u})-J{r',u;). The loss probability of Eq. ^ 
can then be rewritten in the form 


local losses (as in the bulk case) but governed by the self¬ 
interaction process of excitation and back-action. Only 
for certain, rather restrictive simplifi catio ns a viable to¬ 
mography scheme can be formulated®^ the nanoparti¬ 
cles must be small enough such that the quasistatic ap¬ 
proximation can be employed; the plasmonic response 
must be governed by a single plasmonic eigenmode; the 
sinogram must only consist of electron trajectories that 
do not penetrate the particle; the sign of the eigenmode 
potentials must be unique. Although it has been de mon- 
strated that reconstruction is possible in certain cases 
it is obvious that the above restrictions provide a serious 
bottleneck for general plasmon field tomography. 

In this paper we formulate a significantly more gen¬ 
eral scheme, which approaches the reconstruction as an 
inverse problem rather than a tomography scheme. We 
first describe our approach, and discuss possible problems 
and generalizations at the end. First, we decompose the 
dyadic Green tensor into a number of modes Ek(r,uj) 

n 

G{r,r\uj) K.^CkEk{r,u})® Ek{r' ,uj), (6) 
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Pnir,oj) = -^Im|n* • G{r,r,uj) ■ n| 


where Ck controls how much the different modes con- 
J*(r uj) G(r r' w) Jir' uj) \ drdr^^'^^^^^ decomposition. In the following we only 

^ ’ ' ' ’ ’ ’ j consider positions r, r'outside the plasmonic nanoparti¬ 

cle and assume that Ek{r,uj) is a solution of Maxwell’s 
equations. The expansion of Eq. is generally possible 
because G is a symmetric matrix that can be submit¬ 
ted to a singular value decomposition, with Ck being the 
singular values and Ek the orthogonal matrices. In this 
respect, Eq. (§ is similar to a wavefunction expansion 
in quantum mechanics into a complete set of basis func¬ 
tions. 

To be useful as a reconstruction scheme the modes 
Ek {r,u}) should be sufficiently well adapted to the prob¬ 
lem such that a limited number n suffices for a suitable 
representation of G{r,r',uj). Possible modes are quasi 
normal modes of the plasmonic nanoparticleswhich 
have recently received considerable interest, or natural 
oscillation modes of our boundary element method ap¬ 
proach (see Methods). With these modes, the surface 
losses of Eq. become 


( 4 ) 

where dr denotes integration over the spatial variable r. 
Gontrary to Eq. Q , the above expression describes a gen¬ 
uinely non-local self-interaction process where the elec¬ 
tron first induces a field (through excitation of a surface 
plasmon) and then performs work against the induced 
field. 

In Ref. inithe authors tried to interpret Eq. 0 in terms 
of the photonic local density of stated (LDO^ 


( 5 ) 


which is of paramount importance in the field of nanoop¬ 
tics and describes how the decay rate of a quantum emit¬ 
ter located at postion r and with dipole moment oriented 
along n becomes modified in presence of a structured 
dielectric environment. While such interpretation can 
be formally established for nano structures with transla¬ 
tional symmetry along one spatial dimension, it becomes 
problematic for nanoparticles with generic shape.^ 

A different interpretation of Eq. Q in terms of a 
tomography scheme was formulated independently in 
Refs. I8II0I As a preliminary step, let us consider the 
bulk losses of Eq. Q for a given Ri, value. Then, each 
point r inside a medium j contributes with to the 
total loss rate. Within the field of tomographjEIit is well 
known that the three-dimensional profile of 7buik(’') can 
be uniquely reconstructed from a sinogram where bulk 
losses are recorded for all possible propagation directions 
V, using the inverse Radon transform. Such tomogra¬ 
phy reconstruction is significantly more complicated for 
the surface losses of Eq. Q since Tgurf is not the sum of 


hsurf ( Rv ; ^) 
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where {Ra , w) = + vz,uj)dz is 

the averaged mode profile along the electron propagation 
direction. We can now formulate our inverse problem as 
follows. Suppose that one has measured EELS spectra 
Texp for a given loss energy and for various impact pa¬ 
rameters and electron propagation directions. We then 
determine the coefficients Ck such that the entity of mea¬ 
surement data differs as little as possible from the model 
data of Eq. 0, 


. I 

min - 

Ck 2 


^exp {Ri,,Vj) - TsuriiRi,^) 


( 8 ) 
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resulting in a least square optimization (we adopt the 
norm definitions ||a;|||^ = Y,i and \\x\\li = 1^*1)- 

Alternatively, in this work we will use a compressed sens¬ 
ing optimizatioiP^E^ 
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min 

Ck 


rexp(-^'05^) (-^'*>5 


Cfc 
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( 9 ) 

which attempts to minimize the moduli of the expan¬ 
sion coefficients, therefore the scheme is often referred to 
as a Li-optimization, and /r is a parameter that allows 
to switch between genuine compressed sensing and least 
square optimizations.^ For a sufficiently small number 
of expansion modes Ek, the determination of the expan¬ 
sion coefficients Ck is a highly overdetermined problem 
since the measured loss data can be assembled for many 
propagation directions and impact parameters R{, . The 
only pre-knowlege entering our optimization is the self¬ 
interaction-type scattering process of the electron loss, 
Eq. Q, and the assumption that the dynamics of the 
electric fields outside the plasmonic nanoparticles is gov¬ 
erned by Maxwell’s equations. Importantly, once the co¬ 
efficients Ck are determined we have (approximately) re¬ 
constructed the dyadic Green tensor of Eq. Q, which 
allows us to compute all electrodynamic properties in¬ 
cluding the photonic LDOS. 


Results 


To prove the applicability of our reconstruction 
scheme, we generate the “experimental” EELS data Fexp 
using the si mulati on toolbox MNPBEM for plasmonic 
nanoparticles .ISMIl first consider a silver nanorod with 
dimensions 200 x 65 x 30 nm^ and compute the loss spec¬ 
tra for the three selected impact parameters indicated in 
Fig. la. The two prominent loss peaks at low energies 
can be attributed to the dipole and quadrupole plasmon 
modes. Corresponding EELS maps at the resonance fre¬ 
quencies are shown for a few selected electron propaga¬ 
tion directions (rotation angles) in Fig. Ic. The mode 
profiles are reminiscent of the dipole and quadrupole 
surface charge distributions.^^ For the decomposition of 
Eq. ([^ into modes Ek{r,w), we use the information 
about the nanoparticle shape, which in experiment can 
be obtained fro m add itional high-angle annular dark-field 
(HAADF) data,l22llSI and compute the 50 natural oscilla¬ 
tion modes of lowest energy (see Methods). Fig. lb shows 
the modulus of coefficents Ck obtained from either a com¬ 
pressed sensing or least square optimization. Although 
the two approaches give quite different Ck distributions, 
the back-projected EELS maps, obtained by assembling 
the dyadic Green tensor using Eq. ®> and computing 
Fsurf from Eq. both are in almost perfect agreement 
with the original Fexp maps. 

Having obtained the Ck values from the optimizations 
of Eqs. ® we can use Eq. (§ to approximately re¬ 
construct the dyadic Green tensor which allows us to 



FIG. 1: EELS spectra and maps for a silver nanorod. (a) 
EELS spectra recorded at the positions indicated in the inset. 
The peaks at approximately 1.5 eV and 2.7 eV are attributed 
to the dipole and quadrupole plasmon mode, (b) Mode de¬ 
composition of the dipole and quadrupole mode from the col¬ 
lection of rotated EELS maps, using either the least square 
minimization of Eq. ^ or the compressed sensing optimiza¬ 
tion of Eq. §. For each mode the coefficients Ck are normal¬ 
ized to unity, (c) Selected EELS maps for dipole (upper part) 
and quadrupole (lower part) mode and for different electron 
propagation directions (rota tion angles), as computed with 
the MNPBEM toofbox.lS^ (d) Back projected EELS maps 
for the Ck distribution obtained from the compressed sensing 
optimization, using Eq. § for the Green function decompo¬ 
sition and Eq. @ for the calculation of the loss probabilities, 
(e) Same as panel (d) but for Ck distribution obtained from 
the least square optimization. 


compute any electrodynamic response function for the 
plasmonic nanorod. In the following we consider the pro¬ 
jected photonic LDOS of Eq. ([^. Fig. 2 shows the true 
and reconstructed LDOS maps, and compares the qual¬ 
ity of compressed sensing and least square optimizations. 
In particular the inspection of panels (b) and (c), which 
report the LDOS in a plane 20 nm above the nanorod, 
reveals that the compressed sensing results are in very 
good agreement with the true LDOS values, whereas the 
least square optimization completely fails to provide even 
qualitative agreement. This finding seems at first sight 
surprising since both optimization approaches were pre¬ 
viously capable of reconstructing the experimental EELS 
data almost perfectly, as shown Figs. Ic-e. We attribute 
the least square shortcoming to the fact that the EELS 
loss of Eq. 0 is governed by the long-range tails of the 
particle plasmon field distributions, with which the pass¬ 
ing electron predominantly interacts, whereas the LDOS 
of Eq. 0 is governed by the short-range evanescent field 
components. Thus, when the optimization has no strong 
bias on the Ck determination it comes up with the proper 
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FIG. 2: Photonic LDOS of Eq. ([^ and reconstructed LDOS. 
(a) Three-dimensional LDOS distribution, as computed with 
the MNPBEM toolbox (LDOS),^ and the distributions re¬ 
constructed from the compressed sensing (CS) and least 
square (LSQ) optimizations. The projected LDOS pn{r,uj) 
is shown for different projection directions n = x, y, z. (b) 
LDOS density map in a plane 20 nm above the nanoparti¬ 
cle, as reconstructed from the compressed sensing optimiza¬ 
tion. The lower (upper) part of each panel shows the dipole 
(quadrupole) mode, the left (right) part shows the true (re¬ 
constructed) LDOS. (c) Same as panel (b) but for least square 
optimization. The reconstructed least square LDOS has also 
negative contributions, which are set to zero for clarity. 


long-range components, resulting in high-quality EELS 
maps shown in Fig. le, but fails for the short-range com¬ 
ponents which contribute little to the minimization func¬ 
tion of Eq. (|^. In contrast, the compressed sensing op¬ 
timization of Eq. (|^ seeks for a Ck distribution with as 
few non-zero components as possible. For suitable ba¬ 
sis functions Ej. this bias helps to properly select those 
modes which contribute little but still noticeably to the 
loss probability of Eq. Q. We emphasize that such a 
bias for selecting a sparse expansion distribution is by 
no means unique to the problem of our present concern, 
but has been previously highlighted in various studies, 
e.g. in the context of plasmon tomographjffll or single¬ 
pixel cameras,^ and lies at the heart of the compressed 
sensing algorithm. 

An advantage of compressed sensing is that the recon¬ 
struction can in general be performed even with a very 





FIG. 3: Gompressed sensing reconstruction for a strongly re¬ 
duced number of measurement points. The first row shows the 
measurement data for a few rotation angles. In the second row 
we compare the EELS data for a finer sampling mesh (upper 
part of panel) with the reconstructed signal (lower part), find¬ 
ing almost perfect agreement. The last row reports the true 
(upper part of panel) and reconstructed (lower part) LDOS 
maps in a plane 20 nm above the nanorod. 


limited amount of measurement data, and the quality of 
the reconstructed data is usually not strongly affected by 
noise.^In Fig. 3 we show reconstructed EELS and LDOS 
maps for the small number of impact parameters and ro¬ 
tation angles shown in the first row of measurement data. 
As can be seen, the quality of the reconstructed data is 
extremely good despite the limited amount of measure¬ 
ment data. This might be beneficial for EELS experi¬ 
ments which typically suffer from a limited amount of 
rotation angles (missing wedge problem) and where the 
number of measurement points is often kept low to avoid 
sample contamination. 

Finally, in Fig. 4 we compare LDOS maps with recon¬ 
structed maps for a (a,b) bowtie nanoparticle and (c,d) 
cube. For the bowtie geometry we show the LDOS for 
the two plasmon modes of lowest energy, which can be 
labelled as bonding and anti-bonding according to the 
parallel and antiparallel orientation of the dipole mo¬ 
ments of the individual nanotriangles.^ The agreement 
between the true and reconstructed LDOS maps is very 
good, in particular one can clearly observe the strongly 
increased LDOS enhancement in the gap region. For the 
cube we show the dipole and corner modes of lowest en- 
ergyjini finding fair agreement between the true and re¬ 
constructed LDOS maps. We attribute the small differ¬ 
ences to problems of our algorithm when dealing with 
degenerate modes of symmetric particles, which might 
be improved by explicitly accounting for mode symme¬ 
tries.^ 


Summary and discussion 

To summarize, we have shown how to extract the 
dyadic Green tensor of Maxwell’s theory from a collection 
of EELS maps recorded for different electron propagation 
directions (rotation angles). Our reconstruction scheme 
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FIG. 4: (a) True (upper row) and reconstructed (lower row) 
LDOS for a bowtie geometry (total size 215x85x30 nm^ and 
10 nm gap) and for the bonding and anti-bonding modes of 
lowest energy. Color code is identical to Fig. 2. (b) Density 
map of LDOS in a plane 20 nm above the bowtie structure, 
(c) True (left) and reconstructed (right) LDOS for a cube with 
150 nm side length, and for the dipole and corner modes of 
lowest energy.^ (d) Density map of LDOS in a plane 30 nm 
above the cube. 


is based on a singular-value decomposition of the Green 
tensor and a compressed-sensing optimization for the ex¬ 
pansion coefficients. We have demonstrated the applica¬ 
bility of our approach for various elementary nanoparti¬ 
cle shapes. We foresee several improvements for plasmon 
tomography based on EELS. On the experimental side, 
electron holographji^ can provide additional information 
and could allow to disentangle the excitation and mea¬ 
surement channels of plasmonic EELS. On the theoretical 


side, the presented reconstruction scheme works surpris¬ 
ingly well for most nanoparticle geometries, but further 
work is needed to clarify the role of various ingredients. 

First, there are several possibilities for chosing the ba¬ 
sis functions for the decomposition of the dyadic Green 
tensor, Eq. ([^. In this work we have chosen biorthogo- 
nal “constant flux states”l^ that are the eigenstates of 
the Green function evaluated for real frequencies (see 
Methods). They have the advantage that they can be 
computed rather straightforwardly, even in case of de¬ 
generate or near-degenerate modes, on the other hand 
they have to be computed for each loss energy separately 
and several of these modes can govern the plasmonic re¬ 
sponse. Another possibility for a basis are the quasi nor¬ 
mal modes evaluated at the poles of the Green function in 
complex frequency space.EHlll The computation of these 
modes requires an iterative solution scheme,^ however, 
once they are computed they can be used for a large fre¬ 
quency range and in general the plasmonic response is 
only governed by very few of these modes. 

In this work we have considered the situation where 
the basis is already computed for the true nanoparticle 
shape, and have shown that even in this case the EELS 
tomography scheme can be quite tricky. However, our 
approach is less restrictive than it may appear: in princi¬ 
ple, for electron beams not penetrating the nanoparticle 
any basis with modes being solutions of the free-space 
Maxwell’s equations can be employed. Thus, even if a 
slightly different particle shape or dielectric material is 
considered in the computation of the basis, this will not 
necessarily degrade the quality of the reconstruction. In 
this case it might be beneficial to adapt our approach 
such that (i) the modes for the Green function decom¬ 
position are expanded in a given non-ideal basis, and (ii) 
the compressed sensing algorithm seeks for a minimum 
number of decomposition modes. Here it might be ad¬ 
vantageous to use quasi normal modes, because the same 
few modes could be optimized for a whole range of loss 
energies, thus imposing stronger restrictions in compar¬ 
ison to an independent optimization at individual loss 
energies. 

Although further work is needed to establish EELS 
tomography of plasmonic nanoparticles as a robust and 
out-of-the-box scheme, we believe that our present work 
provides an important step forwards for reconstruct¬ 
ing electrodynamic quantities from EELS measurements, 
and makes significant progress with respect to the re¬ 
cently developed tomography schemes that were bound 
to quasistatic approximation and other restrictive as¬ 
sumptions. 


Methods 


Simulations. In our simulation approach we compute 
the LDOS and EELS spectra using the MNPBEM tool- 
bojpSIII] and a silver dielectric function extracted from 
optical experiments.!^ 
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Mode decomposition. For the mode decomposition 
of Eq. (|^ we follow the prescription of Garcia de Abajo et 
al.l^ and compute the natural oscillation modes through 
diagonalization of the E matrix, see Eq. (21) of Ref. HH] 
for details, keeping for the solution of the inverse problem 
the 50 modes of lowest energy. A higher number of modes 
didn’t show a significant improvement in the reconstruc¬ 
tion results. For our mode decomposition it turns out to 
be convenient to use a biorthogonal basis, similarity to 
the quasistatic case.® Our approach closely follows re¬ 
cent related work,® and we introduce the right and left 
eigenmodes Ek{r,uj) and Ek{r',uj) associated with the 
S matrix, respectively. Instead of the decomposition of 
Eq. § we then use 

n 

G{r,r' ,uj) Ek{r,uj) 0 El{r',uj ), 

k^l 

and accordingly also modify Eq. Q. The biorthogonal 
expansion turns out to be advantageous in particular for 
nanoparticles with degenerate modes, as it automatically 
guarantees proper mode orthonormalization. 

Compressed sensing. The least square op¬ 
timization is performed with the built-in Mat- 


lab functions, for the compressed sensing optimiza¬ 
tion we use the YALLl software freely available at 
http://yalll.blogs.rice.edu/. We set the mixing 
parameter /x = 5 x 10“^ and the stopping tolerance 
has a value of 10“^. We take twelve rotated EEL- 
maps for each structure with equidistant angles between 
0 and 180°, each map consisting of 31x51 points. To 
speed up the optimization process we take only 2000 
random measurement points of the generated maps. 
Further, only measurement points with distance more 
than 15nm away from the particle surface are used 
for optimization. For the volume visualization of the 
LDOS we use the MatVTK software freely available at 
http://hdl.handle.net/10380/3076. 
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